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ABSTRACT 



A discrete distribution arises from a Poisson distribution 
with parameter X, when the distribution of X itself is of 
the form cX‘^“^(b-X ) Xe[0,b], where c is a scaling factor 
and a,3,b are strictly positive parameters. However, the 
functional form of the resulting unconditional distribution 
is not particularly tractable, hence the study of the sta- 
tistical properties of the unconditional distribution is 
limited in the study of its mean and its variance. 

In regards to the modelling of a real situation, an 

estimation procedure of the parameters Involved in 
ct““ 1 1 

cX "^(b-X)^ is discussed and a closed form of the proba- 
bility distribution is derived. In addition, when accuracy 
is desired a numerical analysis of the probability distribu- 
tion is also presented. The development of the results is 
continued in Appendix A, as a preparation in computerizing 
the calculation. 

Finally, an application to real data is discussed for 
the purpose of illustrating the model. 
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I. INTRODUCTION 



Let X be a Poisson variate with the probability distri- 
bution depending on a real parameter A, i.e. 

,n -X 

(1.1) Pr{X=nlA=X} = ^^-2- ; n = 0,1,2,... 

n! 

where A Itself is a positive random variable assumed to have 
a prior distribution specified by the probability density 
g(X) over a range of definition R. 



f > 0 V X e R 

g(X) I 

1=0 otherwise. 



/ g(X) dX = 1. 
R 



Depending on the functional form of g(X), the uncondi- 
tional probability distribution of X may result in a wide 
variety of forms. 

With regards to the application the use of (1.1) in the 
field of accident statistics, and acceptance sampling plans 
for manufactured articles based on the counting of defects, 
etc.... is well known. The variations in the accident 
liability from individual to individual in the former case; 
and in the latter an inevitable and continuous change in 
manufacturing conditions, lead to fluctuations in the Poisson 
parameter Involved. 
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However, the true behavior of this random variable Is 
not known. Thus, the distribution of A used so far are 
nothing more than approximations. 

In practice, then. It Is reasonable to estimate the 
first three or four moments of the distribution to get some 
general Idea. 

Referring to accident statistics, this method was ex- 
plained and employed by Newbold In 1927 [!]• The results 
by her data stated that the distribution of A Is 

(1.2) g(i) = 

r(k) 

where k, m > 0. 

That Is a Pearson type III distribution as originally 
suggested by Greenwood In 1919 [1], which Is also known as 
the Gamma distribution with parameters k and k/m. 

It follows that the unconditional probability of X equal 

00 

/ Pr{X=n|A=X}g(X) dX 

0 

(k/m)^ f°°^n+k-l -(k/m +1)X 

J A e dA 

n!f(k) 0 

r(n+k) / k / m 
n! r(k) '^m+k'' ^m+k'^ 



uu n • 



Pr{X=n} = 



(1.3) 
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The random variable X is then Negative binomial distri- 
buted with parameters k and m/(m + k) . 

Using the method of moments, to estimate these parameters, 

/V 

and let k and m are the estimates of k and m; x is the esti- 

2 

mate of the mean, s is the estimate of the variance of X, 
we have 



m = X 



k 



x^/(s^ 



x) 



In regards of modelling a true situation, the above 
model with X given A Poisson distributed and prior distribu- 
tion Gamma (k,k/m) is known as the Negative Binomial Model 
in accident statistics. 

Employing this model, many workers in the accident sta- 
tistics field achieve a good fit to real data. As an example, 
let us consider the following data, referred to accidents 
met at work by 122 shunters during a period of 6 years [1]. 



X 

0 

1 

2 

3 

4 

5 

6 



Observed Frequency 

40 

39 

26 

8 

6 

2 

1 
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We now assume X is Poisson distributed with parameter A 
and Ais distributed Gamma (k,k/m), then apply the Negative 
Binomial model to the data. 

Using the method of moments, with x = 1.27 and s = 1 . 64 , 
to estimate k and m, vie have 

m = 1.27 
k = 4.36 

/s. A 

hence k/m = 3 - 43 > A is then Gamma ( 4 . 36 , 3 * 43 ) distributed. 

Computing the Pr {X=n} by ( 1 . 3 )> and the expected fre- 
quencies for n = 0 , 1,..., 6 results in 



X 


Pr{X=n} 


Expected Frequency 


0 


.3278 


39.99 


1 


.3227 


39.37 


2 


.1952 


23.81 


3 


.0933 


11.38 


4 


.0387 


4.72 


5 


.0146 


1.78 


6 


.0077 


.93 



Performing a chi-squares goodness of fit test, results in 
2 

the value of x as 1.589 with 4 degrees of freedom, hence 
p = S0% . 

However, the Negative Binomial model does not always 
give satisfactory results as in the above example. Since 
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the density function specified by (1.2) is simply an approxi- 
mation to true but unknown situation. In most of the cases, 
the approximations are based on the facts that the random 
variable is positive [2], and that the exponential term is 
included, a simplification is expected in subsequent com.pu- 
tations. Hovjever, the latter factor also makes the assuming 
distribution rigid in the shape. Having tv;o parameters, 
but the shape of the Gamma distribution is mostly dependent 
on one and can only graph three typical curves as illustrated 
in Table 1. 

Thus, (1.2) may not reflect well the true behavior of 
the case where the random variable is distributed in a 
different manner, such as uniformly. 

Working on this assumption, the purpose of this thesis 
is an attempt to look for a more flexible prior distribution 
of the parameter A, then study in detail the resulting 
unconditional distribution of the variate X. 
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TABLE I 



Typical GAMMA curves 
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II. A PROPOSED PRIOR DISTRIBUTION OF X 



Define a density function g(X) of the random variable A 
as follows : 



( 2 . 1 ) 



g(X) 



cx'^'^Cb-x)^"^ ; Ae[0,b] 

^ 0 ; otherwise 



where a>3 and b are strictly positive, and c is a scaling 
factor. g(x) defined as in (2.1), is assumed to have the 
form of the Beta distribution, but the domain of definition 
is [0,b], Instead of [0,1]. For this reason, from now on 
A is said to be Beta distributed. 

A. THE SHAPE OF THE DISTRIBUTION 

Let label (B) be the curve representative of g(X), and 
consider the derivative of g(X) with respect to X. 

( 2 . 2 ) 

(2.3) 

Depending on a and 3, X^^ can be within or outside of [0,b]. 
If X^e[0,b] then g(x) has a maximum (global) or minimum 
(global) in [0,b]. If X^^ [0,b] then g(X) is monotonically 
increasing or decreasing in [0,b]. 



^ g(X) = cx'^ ^(b-X)^"^{(a-l)b - (a+3-2)X>. 
dX ^m ~ ^ a+3-2 
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In general, the shape of (B) near the origin is determined 
by the factor 

Take the limit of (2.1) and (2.2), when X approaches 

zero from the right. It can be seen that (B) will be tangent 

to the vertical axis for a < 1; or intercept this axis for 

6—1 

a ^ 1. By symmetry, the factor (b-X) determines the shape 
of (B) near b. 

Similarly, when X approaches b, (B) v;ill be tangent for 
3 < 1; or Intercept for 3^1 the vertical having the equation 
X = b. 

Thus, for different combinations of a and 3, (B) may 
graph in a wide variety of curves. Some of the typical ones 
are shovm in Table II, Table III, and Table IV. 

B. STATISTICAL PROPERTIES OP THE DISTRIBUTION 

By definition of a density function, g(X) must satisfy 

b 

f g(X) dX = 1. 

0 

or in this case 

c / X°^"^(b-X)^"^ dX = 1. 

0 



Thus, 



r(g+3) 

^a+3-lp(^)P(g) 



and 
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TABLE II 



Typical BETA curves for a < 1 
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TABLE III 



Typical BETA curves for a = 1 
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TABLE IV 



Typical BETA curves for a > 1 
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= IlrteTi 



r(g+B) 



r(a)r(B) 






Let E[A^] be the kth moment about the origin 



E[A^] = 



r(g+B) 



^r(a)r(3) 0 



/ x“-^^-i(b-x)^-^ 



k r(a+3)r(a+k) 
° r(a)r(a+B+k) 



In particular. 



(2.4) 



E[A] = b 



r(a+3)r(a+l) 

r(a)r(a+3+l) 



= b 



a 



a +3 



and 



(2.5) 



pr.2-, _ .2 r(a+6)r(a+2) 



= a(g+l) 

(g+3 ) (g+3+1) 



It follows by (2.4) and by (2.5) that 



( 2 . 6 ) 



Trr.T _ w2 g(g+l) ^2 a 

(g+6)(a+B+l) ‘ 



(g+3) 



of A, 



dX 
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On the other hand, for all teC-tQjt^) the moment genera- 
ting function of the random variable A is 



(2.7) 



M^(t) = 



E[e^^] = 



r(a+3) 



a+6-1 



f(a)r(B) 0 



/ ^(b-X)^ ^e^^ dX, 



which is known to exist but not particularly tractable [3]« 

1 . The Unconditional Distribution of X 

Multiplying Pr{X=n|A=X) by g(X) and integrating over 
[0,b], the unconditional probability of X equal to n is 
given by 



Pr{X=n} = 



r(a+B) 



^r(a)r(6) 0 



b ,a+n-l,, ^ 

/ ^ e'* dX 

n! 



r(g+B) 



n! b”’*'^ ^r(a)r(6) 



I(a,6,n) . 



By letting t = -1, and by substituting a by (a+n) in 
(2.7) and noting that the functional forms of Pr{X=n) and 
(2.7) are now similar. It is then immediate from (2.7) and 
the argument following that Pr{X=n) is not tractable. 

However, the existence of Pr{X=n} can be alternatively 
proved as follows. 

Take the sum of Pr{X=n} for n = 0,1,..., which results in 



Z Pr{X=n} 
n=0 



I 

n=0 



r(g+3) 

bg'^^"ir(g)r(3) 



b 

/ 

0 



x°^~^^-^(b-x)^-^ ^-x 

n! 



dX 



r(g+B) 



/ X (b-X) e E — j- dX 
0 n=0 
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By noting that 



I 

n=0 



n 



n! 



e 



X 



so that 



Z Pr{X=n} 
n=0 



r(g+3) 

b°‘'*'^"lr(a)r(3) 



/ X°^ ^(b-X)^ ^ dX 

0 



( 2 . 8 ) 



= Pr{X£b} = 1 



In addition. 



X°‘‘^^“l(b-X)^ le ^ > 0 Xe[0,b] 

Implying that I(a,3,n) is greater or equal to zero. Then it 
follows by (2.8) that Pr{X=n) exists as a finite positive 
quantity . 

A closed formula of Pr{X=n} will be derived later in 
Section III. 

2 . Estimation Of The Parameters 

Suppose observation is made on a heterogeneous popu- 
lation (S) a mixture of t homogeneous sub-populations during 
a certain period. The observed values of X are recorded, 
and denote Xj^ the maximum of these values. 

Let X^ be the value of X of the ith sub-population, 

and let 
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= max {X^} ; 1 = 1,2,. ..,t. 

and (M) be the corresponding sub-population. 

Clearly, by the expression of g(X), b is identically 
equal to Xj^. However, Xj^ is not knovm so that in practice 
it is estimated as follows. 

Assume for a moment that the members of (M) can be 
Identified from the heterogeneous population (S). Let m 
be the size of (M) and Y be the Poisson variate X, particu- 
larly referring to (M), it is immediate that 



(2.8) E[Y] = Xj^ 



Now let Y^ for i = l,2,...,m be the observed values 
of Y corresponding to the ith individual of (M) . Noting 
that the individuals in (S), with the observed value very 
likely belong to (M) then v;e can assume 



(2.9) ^i “ 

Using the method of moment to estimate X^^, results in 

1 m 

Further, based on the assumption of continuity made 
on the random variable A, we may repartition (S) into t+1 
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homogeneous sub-populations by splitting (M) into two sub- 
populations (Ml) and (M2) in the manner so that (M2) contains 

■f" Vi 

a unique j member of (M) having Estimating A 

in (Ml) and (M2), we have 

^M2 ^ ^ ^M 

A. A. 

Thus, Ajyj 2 > Ajq^ and more general 

^M2 - ^i ^ ~ 1,2, . . . ,M-1,M+1, . . . ,t . 

A 

Hence, Ajy |2 can be used as an estimate of Aj^j, it follows that 
the estimate of b is 

( 2 . 10 ) ^ 

Now let E[X] be the mean, and V[X] be the variance 
of X. By elementary reasoning often used by Bayesians; 

E[X] = E[E(XlA)] 

= E[A] 

V[X] = E[V(XlA)] + V[E(X|A)] 

= E[A] + V[A] 
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Thus, by (2.^1) and (2.5) 



E[X] 



b 



g 

g+3 



and 



V[X] = b‘ 



a(g+l) ^2 a‘ 

- b 



(g-f3,)( g+3+1) 



(g+3) 



+ b 



g 



g+3 



— 2 ^ ^ 

Denote x, s , g and 3 as the estimates of E[X], 
V[X], a and 3 respectively, then 



X = b 



g 



/\ /\ y 

g+3 



and 



( 2 . 11 ) 



g2 ^ -2 g l g+l) _ 



(g+3) (g+3+1) 



- b 



2 g' 



y\ yN o 

(g+3) 



+ b 



g 



g+3 



Subtracting x from (2.11) and dividing by x 



-2 



( 2 . 12 ) 



/N ^ 

(g+1) (g+3) 

/S /N /S 

g(g+3+l) 



-1 




However, 

/\ /s /s 

(2.13) ^ 

g X 

Then, by substituting (2.13) Into (2.18) and adding 1, 
(2.12) becomes 



21 



(2.1^4) 



As As 

b(a+l) 

A 7T. 

x(a+0+l) 



- X -f 
x2 



Now, from (2.13) 



/N /S ^ 1 

J t) 

a + 3 = ot — , 

X 



or 



(2.15) 



As As 



ot + 3 + 1 “ 



a b + X 



Then, by substituting (2.15) into (2.14) to get 



( 2 . 16 ) 



b(g + 1) 

/s 

ba + X 



_ s^ - X + x^ 
^2 



or 



As As A\ As Q 

b(a + 1) x'^ = (ba + x) (s^ - x + x‘^) 



Solving for a in the above equation, it comes out 



(2.17) 



x(bx + X - x^ - s^) 
a = K 5 ;;; 

b(s‘^ - x) 



Substituting the value of a into (2.11), to have 



( 2 . 18 ) 



B = a(^ - 1) 

X 
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Now observe that estimating by by implies 

(2.19) b > X or > 1. 

X 

Furthermore, considering (2.9) > this equation is symbolically 
written as 



V[X] = V[A] + E[X] 

Noting that the variance of a random variable is always 
positive, then 



V[X] - E[X] = V[A] ^ 0, or equivalently to 



( 2 . 20 ) 



2 

s - X 



= b 



2 a(a + 1) 

— 7^ — 7^ 7N — r: 

(a+B) (a+3+l) 



- b‘ 



a 



(a+3) 



2 - 



> 0 



Also observe that by (2.17), (2.l8) and (2,19), if 
a is negative (positive) it follows that 3 is negative 

/V /V 

(positive) and vice-versa. Suppose now that a and 3 both 
be negative. This implies that 



a + 1 



a + 3 + 1 



< 



/s 

a 

Tv 7T 

a + 3 



since 



/N 




a+3 
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and leads to 



/\ /N /\ Q 

a(a +1) a 

(a + 3)(ot +3 + 1) (a + 3) 

As a result, (2.20) is negative, which is a contradiction, 

/s /s 

and thus a and 3 should both be positive. 

/N /S 

Now assume a is equal to zero, 3 then is also equal 

to zero by (2.l8). This leads to either or both of the 
following cases. 

a. X = 0 

This is impossible, since the positive random 
variable X can not have zero mean. 

b . X = b 

Which contradicts to (2.19) 

A A 

Thus, a and 3 should be both strictly positive as 



required in (2.1). 



III. CLOSED FORMULA OF PR{X=n} 



As mentioned in the last section, the exact mathematical 
expression of the Pr{X=n} can not be derived by the present 
method of integration. Approximation is then used to esti- 
mate through a closed form of the integral I(a,3jn). 

First, by noting that 



e 




k=0 . k! 



Hence , 



(3.1) 



Ka, 3,n) 







Now, consider the Taylor’s expansion of 




(3.2) 




+ E (-1) 
k=2 




Multiplying (3.2) by f(3)b 




(3.3) 
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Then, suppose for some values of a and B, the following 



equation is known to approximately hold 



(a+B+n)^ = ( a+B+n) . . . ( a+B+n+k-1 ) 



it follov;s then that 



(3.4) 
and thus 

(3.5) 



r(a+B+n) (a+B+n)^ 
( 3 . 3 ) becomes 

r(B)b“+2-'^-i z 

k=0 



= r(a+B+n+k) 



/_,Nk b^r(g+n+k) 

^ ^ k!'r(a+B+n+k) 



Taking into account (3.4), then comparing (3.5) with ( 3 . 1)5 
it leads to 



(1 + 



b x-a-n 



a+B+n 



) 



[r(B)b 



a+B+n -1 r(a+n) 



-1 






recalling from the last section, Pr{X=n} is written as 



Pr{X=n} 



r(g+B) 

n!bCt+3-lr(a)r(B) 



I(g,B5n) 



Hence , 



(3.6) Pr{X=n} 



b^r(a+B)r(g+n) /, b x -a-n 
n!r(a)r(g+B+n) ^ g+B+n^^ 
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It can be seen that formula (3-6) above is an overestimate 



of Pr{X=n}, i.e 



\ _ b^r(a+B) r(a+n) ^ b ^-a-n 

- nirralT'Ca+T+n)' ■ 



a+B+n' 



- e 



n 



where is a positive quantity, depending on a,B and n. 
In fact, let 



(3.7) T 



k 



b^r ( g+B) b^r ( g+n+k) 



1 nTTTa) i^ir(ot+g+n)(a+B+n)^ 



(3.8) 



nk _ b^r(g+B) b^r( g+n+k) 



n!f(g) k!r(g+B-'-n+k) 



Multiplying both numerator and denominator of by 
(g+B+n) . . . (g+B+n+k-1) it results in 



(3.9) 



,k _ b^r (g+B) b^r(g+B+n) (g+B+n) . . . (g+B+n+k-1) 



n!f(g) 



k!f(g+B+n+k) (g+B+n) 



k 



Taking the difference 5^^ of (3-9) and (3.8), and noting 
that for 6g = 6^ = 0 then for k ^ 2 



. _ b^r( g+B) b^r( g+n+k) r (g+B+n) . . . (g+B+n+k-1) , ■. 

*^k n!f(g) k ! f ( g+B+n+k) ^ (g+B+n)^ ~ 



k. 



= b^r(g+B) b^r( g+n+k) . j x _ . 

n!f(g) k!f(g+B+n+k) g+B+n'' 
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Thus, for k ^ 2, 6 ^ increases as long as the numerator 
Is larger than the denominator, then for some k, it converges 
to zero. The speed of the convergence is proportional to 

(X + 3 ♦ 

It follows that 



n 



00 

z 

k=2 



(- 1 )^ 6 , 



Hence, for a + 3 large enough, is small, (3-6) is then a 
good approximation of Pr{X=n}. 

In practice, (3*6) can be applied vastly together with 
some correction. Before doing so, however 6^ should be 
expressed more explicitly. 

Consider 



k-1 

n 

j=0 



(1 + 



a+3+n 



1 , 



and assume that — is small enough so that 

a+3+n 



k-1 

n 

j=0 



(1 + 



a+3+n" 



1 



a+3+n 



2 

a+3+n 



+ 



. . . + 



k-1 

a+3+n 



(k-l)k 

2(a+3+n) 



hence , 



. ^ b’^r(a+3)r(g+n) 

k 2(a+3+n)n! r(a)r(a+3+n) 



/b^(a+n) . . . (a+n+k-1) 

^k! (a+3+n) . . . (a+3+n+k-l) 



(k-l)k) 
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and 



e 



n 



b’^r(a+B)r(g+n) 

2 ( a+6+n)n ! r ( a) F ( a+6+n) 



00 

1 (- 1 )^ 
k=2 



b^(g+n) . . . 
k ! ( g+6+n) . 



(k-l)k } 



The factor in brackets does not depend greatly on n. 
approximately has the common value a for some small 



P = b^r(g+B)r(g+n) ^ 

n 2(g+B+n)n! f(g) F (g+B+n) 



now suppose there is an N, such that 



S = 



V b'^r( g+B) r ( a+n) 
^^On!r(a)f(g+B+n) 



(1 + 



g+B+n 



-g-n 



b. Pr{X=n} is close to zero for all n > 



Then, 



N 

S - 1 = Z e 

A n 



^ ^ b^r(g+B)f(g+n) 

^_Q 2( g+B+n)n ! f ( g) f( g+B+n) 



Solving for a in this equation, it results in 



(g+n+k-1) 

. . (g+B+n+k-lO 



It 

n . Thus , 



1 
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N 

a = (S-1) [ E 
n=0 



b^V(a+Q)T(a+n) 

2 ( a+6+n)n ! F ( a) F (a+6+n) ■' 



Substituting o back into (3»10), to get the approximative 
e^, and finally, with the correction, (3.6) is modified as 



(3.11) Pr{X=n} 



b^r(a+3)r(g+£n) r/_ , b s-a-n 

n! r(a) r(a+3+n) ^ a+g+n"^ 



g 

2(a+3+n) 



} 



The application of the above results will be seen in 
Section V. As is shown, the estimate of Pr{X=n}, using (3.11) 
agrees to tv/o decimal digits, in overall, with the results 
using (1.3)j in Section I, for b = 6 and a + 3 = 13.009782. 
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IV. NUMERICAL TECHNIQUE IN APPROXIMATING PR X=n 



In the last section, a closed formula of Pr{X=n} was 
derived. Although it proved itself a fairly good estimation, 
the error Involved was uncontrollable, and thus the validity 
of the proposed distribution of A in Section II cannot be 
judged. Hence, numerical techniques of integration are 
suggested to evaluate the Pr{X=n} through the integral 
I(a,3>n). Several methods have been considered, such as 
Gaussian quadratures and Romberg's extrapolation to the 
limit [4]. However, the direct application of these methods 
leads to unsatisfactory programming considerations and 
unbounded error estimates. 

Remember that Pr{X=n} can be written as 



Pr{X=n} 



b^r(g+B) 

n! r(a) 



y ( b^r(a+n+k) 

^ k!r(a+B+n+k) 



equivalent to a series of alternate sign terms which converges 
monotonically to zero for large enough k. This behavior 
guarantees that Pr{X=n} can be estimated by summing the 
first K terms, and the absolute error is no larger than 

S t/ 

the absolute value of the K + 1 term, l.e. 

b^r(g+B)b^'*'^r(a+n+K+l) 

n!r(g)(K+l) !r(g+B+n+K+l) ‘ 
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It seemed to be a powerful method, however the conver- 
gence is rather slow. Moreover, if the estimation is done 
by means of a computer, the precision will be altered by 
the propagation of the error in repeated subtraction of 
nearly equal terms in the alternating series. 

Thus, a special technique is required v/hlch is expected 
to Improve the convergence and also eliminate the error 
propagation. 

This technique is simply a variant of the Gaussian 
quadrature method. 

A. CONVERGENCE ACCELERATION 

_X 

The idea is to approximate the factor e in the integral 

Pr{X=n} = / x“'*'^"^(b-X)^"^e"^dX 

n!b“ P -^r(a)r(3) 0 

with an orthogonal polynomial p (X) of degree s by least- 

s 

squares approximation method [4]. In other words, it is to 

find p (X) so that 
^s 

-X 2 

(4.1) E(s) = / {e - p (X)} dX is minimized. 

0 ® 

From the theory of this method p„(X) can be expressed 
as 

s 

(4.2) p^(X) = Z d,P,(X) 

s j=0 J 
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where the {P.} is a sequence of ortho-polynomials of degree 
J 

j and satisfies the properties below; 



(^.3) 






PqCX) = 1. 



b (S ; if k = j. 

/ ?AX) P.(X) dX = { J 
0 ^ 'J lOjifk/^j. 



d ,= ^ / e ^ P,(X) dX 

J Sj 0 J 



P 



J+1 



This method also provides a recurrence formula of 

(X) in terms of P.(X) and P._^(X) which is 

J J -L 



C^l.5) 




(X - B ) P (X) - CP (X) ; 

J J J J 



where B. and C. are defined as 
J J 

1 ^ 2 

(^1.6) B. = ^ / XP.(X)^ dX 

J 0 J 



and 



(4.7) 




for j > 1 

for J = 0 



For the sake of simplicity, let now denote P.(X) by P. and 

J J 

P„(^) by simply p . Then, taking into account (4.1) and 
s s 

(4.2), E(s) was expressed as 
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^ -X 2 

E(s) = / {e - p } dX 

0 ^ 



-2X 



s s 



= / e dX + E E d.d / p p dX - 

0 j=0 k=0 J ^ 0 ^ ^ 



2 L d, r -Xp ,, 

J=0 ® 



by (^*3) and E(s) is reduced to 



E(s) = 1/2 (1 - e"^^) - E d.^ S. 

j=0 J J 



Thus, for a value of b, E(s) is decreasing when p takes 

s 

a higher degree, and as b increases s should be also increased 
in order to keep E(s) as small as desired. 

As an example, some of the values of E(s) v/ere displayed 
in Table V. 



TABLE V 



The least-squares error E(s), for b = 6 



E( 0 ) 




.33415549 




E(l) 


= 


. 109724^17 




E( 2 ) 




.20317197 


10 ^ 


E( 3 ) 


= 


.2362^1092 


10'2 


E( 4 ) 


= 


.18656389 


10"^ 


E( 5 ) 


= 


.10601508 


10"^ 


E(6) 


= 


.^ 152841^17 


10“^ 


E( 7 ) 


= 


.15041719 


lo""^ 


E(8) 


=r 


.39915586 


10"^ 


E( 9 ) 


= 


.86497233 


H ( 
H . 

1 

0 

t — 1 


E( 10 ) 


= 


.15588653 


10 
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B. GENERATION OF THE SEQUENCE {P^} 

By definition, can be written as 

(4.8) P = Z ; 

J k=0 ^ 

1 k 

where the ‘ s the constant coefficients of the X *s. 

It is now to generate explicitly the i.e., to find the 

Aj^*^ • s in terms of b. Before doing so, let us express , 

B . , and d. in a computable form. 

J J 

Taking into account (4.8), (4.3) is then written as 

(4.9) S = E E / A ^A dX 

J k=0 £=00 ^ ^ 



J 


j 


A 




b 




= E 


E 


k 


* 






k=0 


z=o 




k+£+l 


0 




j 


j 


A, J 


^ j^k+£+l 






= E 


E 


k 


it 






k=0 


o 

II 




k+£+l 






Similarly, (4.6) becomes 








(4.10) B. = :^ 


J 

E 

k=0 


j 

E 

£=0 




1 


j 


j 


A j, J.k+t+2 


b 


Q 


E 


E 


k £ 






J 


o 

II 


II 

o 


k+it+2 




0 




J 


i ^ 


J A 0vk+£+2 




Q 


E 


E 


k it 


• 






o 

II 


II 

o 

1 


k+£+2 
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and finally (4.4) becomes 



1 ^ 

(4.11) d = ^ E A, 

J k=0 ^ 



In particular, for j = 0 applying the results just 
developed to get 

Sq = b ; Bq = b/2 ; and d^ = (1 - e"'^)/b 
Then, by (4.5) 

= (X - b/2) . 

With the expression of as above, applying again (4.9), 
(4.10) and (4.11) for j = 1, it results in 

= h^/12 ; = b/2 and 

d^ = -e"^(12/b^ - 6/b^) + 12/b^ - 6/b^ 



r -X 
Jab dX 



/ ^ -b k ! .111 

‘-/o ® Tk^ ^ 



For j ^ 1, P-i-i is generated recursively by (4.5) and 
by means of a Fortran-Pormac program [5]> vfhich has the 
capability of manipulating symbolic mathematical expressions. 
This program also computes and directly by 
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(^.9)) (^.10) and by (iKll). Some of the first , so 
generated, are displayed in Table VI. 

The following is the Algorithm, v/hich the Portran- 
Formac program mentioned above is based on. 

ALGORITHM 1. 

Step 1. Read. in P- , S^ , , d^ 

0 0 0 0 

Read in P^, S^, C^, d^ 

Set j = 1 

Step 2. Compute P-.t by (4.5) 

Compute by (4.9) 

Compute (4.10) 

Compute + by (4.11) 

Step 3. Updating j by 1 
Go to Step 2. 

With Portran-Pormac programming, given P*_-i 5 Pj > we can 
obtain Pj+q_ et any desired j. However, in practice, if b 
is known based on (4.5), the coefficients Aj^*^ ' s can be 
obtained with Fortran programming alone as follows. 



= -B A^' 



c aJ"1 • 

j 0 



for 1 < m < j-2 



A'^'*'^ = -B.A^* - C.A'^ ^ + A^ , ; 

m J m J m m-1 



and 



= 

J-1 

A^-^l = A 



- B.A-! T + A^. ^ 
3 J-1 J-2 






j j-1 * 
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I 



TABLE VI 

Pj and dj for j = 2,3. 

P2 = -b^ +b^/6. 

P^ = -3bX^/2 + 3b^X/5 - bV20. 

d^ = -e"^(360/b^ + l80/b^ + 30/b^) + 360/b^ 

- 180/b^ + 30/b^. 

d^ = -e"^(l6800/b'^ + 8400 /b^ + l680/b^ + l 40 /b^) 

+ 16800^"^ - 8400 /b^ + l680/b^ - l 40 /b^. 
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But, the A • s computed with Fortran programming v/ere 

not as precise as those derived by Fortran-Formac , since it 

was dealing with a large amount of computation, while the 

• s are evaluated directly from symbolic expressions in 

Fortran Formac programming. However, the differences v/ere 

negligable when j was less than 10. 

Thus, taking account of (4.8) and with the Aj^*^ ' s found as 

above, p is now written as 
s 

2 <3 • v 

(4.12) Pc = 2 2 

^ j=0 J k=0 ^ 



Now, let P(n,s) be the estimate of Pr{X=n}, approximate 
e~^ with p . Taking account of (4.12), 

o 



P(n,s) 



E d, E A^. 

n!b“ r(a)r(3) j =0 ‘^k =0 ^ 



0 



Further, 



let t^' 
k 



d. A^ 
J k 



Then, 



P(n,s) 



b^r(g+B) ? i .j ^k r(g+n+k) 

n ! r(g) j^_Q ^k f(g+B+n+k) 



and P(n,s) above was computed by means of an iterative 
method, based on the following algorithm. 



ALGORITHM 2. 



Step 1. 



Set J = 0 
Let P(n,s) 



_ b’^r( g+g) ^o r(g+n+k) 
n ! r (g) ^0 r ( g+3+n+k) 
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step 2. Set j = J + 1 



Compute T( j ) 



^ b^r(g+e) ^ .3 , k r(g+n+k) 

n ! r ( g ) r(g+B+n+k) " 



Set P(n,s) = P(n,s) + T( j ) 



Step 3. If |TCJ)| > r Go to Step 2 
If |T(J) I < r Stop 

It is noted here, that Step 2 was called the J + 1^^ 
iteration. The stopping criteria in Step 3 is only for the 
purpose of fixing the idea, it was replaced later by a more 
appropriate one discussed in the next subsection, when the 
error in P(n,s) is discussed. 



C. ERROR ESTIMATION 

Let e(n,s) be the error committed in approximating P(n,s) 
to Pr{X=n}. Thus, 



(il.l3) e(n,s) = Pr{X=n} - P(n,s) 



r(g+B) 



n!ba+B-lr(a)r(B) o' 



,g+n-l,. .sB-l, -X ^ 

X (b-X)^ (e -Pg)dX 



Now by the mean value theorem of integration, and by 

noting that X ~ (b-X)^ does not change sign in [0,b], 

then there exists a X*e[0 ,b] , which depends on X°^’'’^”^(b-X) ^ ^ 

and on e”^-p (X), so that 
s 



(4.l4) e(n,s) = 



{e“^*-p ( X*) 

o 



r(g+B) 



n!b 



a+6-1 



r(a)r(a)o 



/ x“'^’^"^(b-x)^ ^dX 



In general, such a X* in (4.l4) is not known, hence e(n,s) 
cannot be found. However, by using the Cauchy-Schvjarz 
inequality the bound u(n,s) of the error is derived as follows. 



I e(n, s) I £ u(n, s) 



r(g+B) 

njba+3-lr(a)r(g) 



0 



[ /(e'^-p dX]2 
0 ® 



Recalling that. 



E(s) 



-X v2 ,, 1 ,, ^-2bv 

/ (e -p ) dX = 5- (1-e ) 

0 ^ ^ 




and that 



0 



. 2(a+6+n-l)-l f ( 2g+2n-l) f ( 2B-1) 
r(2g+ 2S+2n-2) 



then 



C4.15) uCn,s) 



b”"2r(g+B) rr(2g+2n-l)r(2B-l) r-/„^-l2 

n!fCg)r(B) F ( 2g+2B+2n-l ) 



On the other hand, the error e(n,s) in (4.13) can be 



alternatively written as: 



I 



e(n,s) 



r(g+g) 

n!^a+B-lr(^)r(3) 



0 



r(g+B) 

njba+B-ir(g)r(e) 



. -.g+n-1/, -inB-I 
/ X (b-X) p dX. 

0 ® 



Now, applying the mean value theorem of Integration to the 
two above integrals, and let 



V 



r(g+B) 

njba+B-lr(g)r(B) 



X“'^^"^(b-X)^"^ dX 

0 



b”r(g+B)r(g+n) 
n!r(g)l'(g+B+n) ’ 



then 



(iJ.l6) e(n,s) = {e ° - p^(X^)} V 



where Xp. and X, e[0,b]. 

-h 

Adding and subtracting to (from) (4.l6) the quantity e V, 
we have. 



(4.17) e(n,s) 




(X^)} V + {e 





} V 



-X. 



= e 



V - P(n,s) + {e 



-X. 




} V 



Further, assume that X^ ~ then (4.17) becomes 
(4.18) e(n,s) = e ^ Y - P(n,s) . 
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It can be seen that (4.l8) Is a good estimate for e(n,s) 



only when the above assumption holds. 

The in (4.l8) can be found by solving 

(4.19) f(X) = Pg(A) - ^ Q 

Based on the above analysis, algorithm 2 is now modified 



as , 




Step 1. 


Set j =0 

p C n " ^ ~ ^ f(oi+3) ^0 r((x+n+k) 

t n!r(a) 0 r(a+6+n+k) 


Step 2. 


Set j = 0 * +1 

rr,(.s _ b^r(a+6) 1 r(a+n+k) 

Compute T(j) n!r(a) r(a+B+n+k) 

P(n,s) = P(n,s) + T( j ) 


Step 3- 


Compute u(n,s) 


Step 4. 


If u(n,s) £ r Stop 

If u(n,s) > r Go to Step 2 



However, for the purpose of comparison, the approximate 
e(n,s) as in (4.l8) is also computed, and it is suggested 
that (4.19) be solved by Newton iterative method which is 
suitable for this problem [4]. In fact, f(X) is continuously 
differentiable in [0,b] as the Newton method required. In 
addition, the derivative of f(X) is easily computable. 
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Taking into account (4.12), then 





f(X) - E E t.-^ X^ - — 

j=0 k=0 I' '' 


it follows 


that , 


(4.20) 


f(A) = E E k tj 

J=1 k=l I' 



Moreover, (4.20) evaluated at does not vanish as p (X) 



approaches 


-X 

e 


On the 


other hand, a starting point X^*^^ for the Newton 


iterations 


is nicely provided by: 


(4.21) 


-X^°^ _ P(n,s) 
e - V 



In such conditions, the Newton methods are very efficient, 
i.e., converge rather fast (quadratically) to X^. 

For a detailed treatment of the analysis, the reader is 
referred to Appendix A. 
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V. AN ILLUSTRATION 



Consider the following example. 

A sample of 621 white male children is observed, when 
they were from 8 to 11 years old. The data reported here 
is a part of a large study of childhood accidents conducted 
by the State of California Department of Public Health. 
Subjects to study were those children v;hose families sub- 
scribed to the Kaiser Foundation Health Plan [6]. 

The number of injuries and the observed frequency were 
tabulated below. 



X Observed Frequency 



0 


2^10 


1 


192 


2 


107 


3 


52 


4 


17 


5 


9 


6 and more 


4 



For purposes of comparison, let us first apply the 
Negative Binomial to these data, i.e. assuming that X is a 
Poisson variate with parameter A, and A Itself is Gamma 
distributed . 
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A. THE NEGATIVE BINOMIAL MODEL 



1 . Distribution of A 

Using the standard method to estimate the mean and 
the variance of X results in 

X = 1.13 

and 

s^ = 1.52 

Hence, by (1.^) and (1.5), 
m = 1.13 



and 



k = 3.24, thus k/m =2.8? 

It follows that from (1.2), g(X) is Gamma (3.24,2.87). 
The mean and the mode are respectively, 

E[A] = 1.13 

X = .785 

m 

The curve representative (C) of g(X) was shown in 



Table VII. 



2 . The Unconditional Probability of X = n 



With 


the values of 


k and as above, and 


Pr{X=n> is 


r(3.24 + n) 


.15 X .26“^ 




n! 


Then, for n = 


: 0,1, 2, 3, 4, 5, 6 


and more : 


n 


Pr{X=n} 


Expected Fr^ 


0 


. 38063 


236.37 


1 


.31375 


197.38 


2 


.17374 


107.88 


3 


.07836 


48.59 


4 


.03149 


19.55 


5 


.01176 


7.30 


> 6 


.00627 


3.89 


2 

X = 


1.179 with 4 


degrees of freedom 



75% < p < 90% 



B. APPLICATION OP SECTIONS II, III, IV. 

1 . Distribution of A 

Assume now g(X) Is of the form (2.1), and assume 
further that the probability of X is greater than 6 is small, 
then by (2.9) 
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b = 6 



(B) 



tcj 



cr; 

.o 



C-J 



CJ 

C3 




HI Z 
C-J c_ 



CO uO 



'21 21 
“D ID 



o 



0^1 — I lO 
I I 

Li-- Li-- ^ ■ — 3 

- 'r cr 

“ ^ LLJ 

.■'■ .’■ I . ; — % 

LJ-- Ui-- *— I 1 / 



Cl CT 



; r>^-> 

! I I • 






l± 



1 I 
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It follows by ( 2 . 17 ) and by ( 2 . 18 ) that 



a = 2 .i|il 0642 

B = 10.569140 

and, 

a + B = 13.009782 



hence , 



g(X) =..13 10 -« ( 6 -x) 5 - 5691 '< 



with mean and mode respectively. 



E[A] = 1.13 



X = .785 

m 



The curve representative (B) is shovm in Table VII. It 

is noted that (B) and (C) have the same mode, X = .785. 

m 

2 . The Unconditional Probability 
a. Using Closed Formula. 

Substituting the values of a and B found above in 

(3.6) 



- 6^ 3.85 10^ r(2. 440642 + n) 

^ n! 1(13.009782 + n) 



(1 + 



■) 



- 2 . 440642 -n 



13.009782 + n 
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Then, for n = 0, 1, 2, 3> ^ 
of the probabilities of the 
are the following 

n 

0 

1 

2 

3 

4 

5 

6 



5 and 6, the rough estimates 
unconditional of X equal to n 

Pr{X=n} 

.39629 

.33017 

.18628 

.18628 

.08685 

.01338 

.00463 



Now applying the correction to the above results. 
By (3.11) and (3.12), 

s = 1.0533 
s - 1 = .0533 

and 

a = .4001 



Thus, the results should be corrected as 
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n Corrected 

0 .38091 

1 .31049 

2 .17522 

3 .08071 

4 .03279 

5 .01211 

6 .00414 

It can be seen that with the correction, i.e., 
using formula (3»13)» the estimates of the probabilities 
of X equal to n, can be compared within two decimal digits 
to the results found by the Negative Binomial model. However, 
the significance of these figures is uncontrollable, hence 
there is no use to perform a goodness of fit test here, 
b. Using Numerical Technique 





Fixing in advance r = .156 10 , 


then after 11 


Iterations 


P(n,s) . converge 


to the following 


figures . 


n 


P(n,s) 


u(n,s) 


Expected Frequency 


0 


.383299 


2.0 lO"”^ 


238.03 


1 


.313539 


2.9 lo"”^ 


194.71 


2 


.173662 


2.1 lO"”^ 


107.84 


3 


.079539 


1 

0 

1 — 1 

OJ 

• 

1 1 


49.39 


4 


.032189 


6.2 10“® 


19.39 


5 


.011869 


2.8 10"® 


7.37 


6 


.004057 


1.1 10"® 


2.51 
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Based on the above results, and by noting that 
the sum of P(n,s) from n = 0,1,... equals to one, a chi- 
square goodness of fit test is performed, the results are 
= 1.035 with 3 degrees of freedom and 75^ < P < 90 %. 

It is noted that the values of the P(n,s) and 
u(n,s) are rounded off to the significant digits. For the 
convergence and the actual values of P(n,s), u(n,s) and 
e(n,s) the reader is referred to Appendix B. 

The summary of the results obtained in this 
section are shown in Table VII. 
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TABLE VIII 



SUMMARY OF RESULTS 



X 


(i) 


(ii) 


(ill) 


(iv) 


(v) 


(vi) 


( vll) 


0 


. 38063 


. 38091 


.3832990 


240 


236.37 


236.54 


238.03 


1 


.31785 


.31409 


.3135385 


192 


197.38 


195.04 


194.71 


2 


.17374 


.17522 


.1736619 


107 


107.88 


108.81 


107.84 


3 


.07826 


.08071 


.0795389 


52 


48.59 


50.12 


49.39 


4 


.03149 


.03279 


.0321885 


17 


19.55 


20.36 


19.98 


5 


.01176 


.01211 


.0118698 


9 


7.30 


7.52 


7.37 


6 


.00627 


.00417 


.0059066 


4 


3.89 


2.61 


3.66 



Note : 
( 1 ) 
(il) 
(iii) 
(Iv) 
(v) 

(Vi) 

(vii) 



Pr{X=n} (negative binomial) 

Pr{X=n} (closed formula 2 with correction) 

Pr{ X=n} (numerical technique) 

Observed frequency. 

Expected frequency (negative binomial) 

Expected frequency (closed formula 2 with correction) 
Expected frequency (numerical technique) 
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I 

I 

I 
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VI. CONCLUSIONS 



As is shown in Section V, an application to the real 

situation of the new model is presented. In it we assume 

that the data is Poisson distributed with parameter A and 

the parameter itself is also a random variable with density 

function of the form as specified in (2.1). The validity 

of the model is then made by performing a chi-squares good- 

2 

ness of fit test. As the results, x = 1*035 with 3 degrees 

of freedom and 75% < P < 90%. 

On the other hand, modelling the set of data as Negative 

binomial model, and using the same validation procedure the 
2 

results are x = 1*179 with 4 degrees of freedom and again 
75% < P < 90%. 

Comparing the two results, and referring particularly 
to the expected frequencies, it can be seen that those of 
the new model are closer to the real data. However, because 
the distribution specified in (2.1) has 3 parameters, i.e., 

1 parameter more than the Gamma, causing 1 degree of freedom 
less in the chi-squares goodness of fit test the results do 
not differ significantly. Availability of an extra parameter 
is one of the major advantages of the new model in regards 
to the modelling of any real situation. However, using the 
maximum of the observed values of X in a single period as 
an estimate of the parameter b in Section II, is not an ade- 
quate method. It leads to difficulties in extending the 



5 ^ 



model in order to predict the distribution of X in the 
next period. 

Moreover, the exact mathematical expression of Pr{X=n} 
cannot be derived, hence the knowledge about the statistical 
properties of X is limited to the mean and the variance. 

In Section III, a closed formula of Pr{X=n} is then 
derived. However, it can be applied only in cases where 
the sum of a and 3 is large enough. 

Finally, a numerical method for estimating Pr{X=n} is 
presented in Section IV. It can be seen that this is a 
powerful method when high accuracy is desired. However, 
it is also complicated and can be best used by means of a 
high speed digital computer. 

As recommendations for further investigations, the 
following are considered as potential topics: 

1. Determining the exact expression of Pr{X=n}. It 

is believed that such a formula can be obtained by advanced 
mathematical analysis. 

2. Determining an adequate estimation procedure for 
the parameters Involved in the prior distribution of A as 
specified in (2.1). 

3. Continuing with an extension of this model to a 
bivariate or more general, k-varlate model. 
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APPENDIX A 



COMPUTERIZING THE CALCULATION OF P(n,s) AND 
DETAILED TREATMENT OP THE ERROR ESTIMATION 

The following shows the development of the modified 
algorithm 2 in Section IV into details to computerize the 
calculation of the P(n,s) and the error involved u(n,s) as 
well as the estimated error e(n,s). 

Suppose the values of a,3jb and Pg, dg, Sg, Bg and P^, 
d^, S^, were already input. We now wish to compute P(n,s), 
u(n,s) and the corresponding e(n,s) for n = 0,1,..., N. 

Since the computation is done by an iterative method, the 
stopping criteria should be determined first. 

Let r* be the maximum of the u(n,s) for n = 0.1..., N 
at the s iteration. The iteration is then stopped when 
r* < r, where r is a small quantity determined in advance. 

Step 1. Initialize all the P(n,s) to zero. 

Set E(s) = .5(l-e"^^) 

Set j = 0, where j is the subscript of the P.. 

J 

Step 2. For j = 2,3j«** and note that s is the actual 

value of J in this step, then 

Compute P., S , B , d by. (4. 5), (^-9), (^.10), 

J J J J 

and (4.11) 

Compute t^'^ . 

Set E(s) = E(s) - (d.)^ S. 

J 0 
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step 3» For n = O,!,..., N and for j = 0,1,... 



Compute T(j) 



^ b^r(g+g) i j k r(g+n+k) 
n ! r( g) j^_Q ^k T ( g+B+n+k) 



(a.l) Set P(n,s) = P(n,s) + T(j) 

Compute u(n,s) by (4.15). 

Compute e(n,s) by (4.l8) for j > 0 
First by solving for in: 

(a. 2) f(X) = Pg(X) - = 0 



As mentioned previously (a. 2) is solved using Newton 
iterative method [4]. The idea is to generate the sequence 
where i = 0,1,... by the iterative formula: 

(a. 3) = a'-") - 

starting from X^^? as defined in (4.21). 

Furthermore, if the sequence so generated converges to 
some point l.e., for some inter L, we have 

(a.i.) = J 



so that (a. 3) now becomes 



then it is immediate that f(?) = 0, and that C is the 
desired 

Preliminary, by (a.3)> f(X) should be continuously 
differentiable and its derivative with respect to X should 
not vanish at ? = X^. However, referring to the analysis 
done in Section IV, it can be seen that fC^) satisfy both 
of the above conditions. 

On the other hand, starting with X^*^^ as found in (4.21) 
the sequence will converge to X^ quadrat Ically . The 

proof is as follows. 

Consider, 

(a. 5) h(X) = X - 



Taking the derivative of (a. 5) with respect to X, 



(a. 6) 




(Pg(X) - f"(X) 

[f(X)]2 



Now, by noting that, in the vicinity of X^, (a. 6) is 
small with respect to 1, and it vanishes at X^, then making 
use of the discussion in page 57 of the reference, it follows 
that the iteration converges quadratlcally to X^. 

But, the integer L as in (a. 4) is not known, and usually 
is estimated by L* in advance. Hence, X^ may be found 
exactly or with some tolerance within L* Newton Iterations, 
if X^^^ is close enough to X^. 
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Based on the above argument, we may then stop the 
searching process for if for some m the following conditions 
are met. 



(a. 7) 



(a. 8) 



|^(m+l) _ ^(m)| ^ ^ 

I < e 



It is noted that m £ L* and that <}), 0 are small positive 
quantities . 

The above conditions are checked by associating them 
with the binary function Flag defined as below 

^0 if (a. 7) or (a. 8) are satisfied 
Flag = \ 

L 1 otherwise 



Then during the course of searching for a X^, if Flag 
comes out with the value 1, either one of the following 
corrections have to be made. 

1. Increasing L* 

( 0 ^ 

2. Substituting X^ by X^ , 

The corrections may be repeated until Flag = 0. 

Summarizing the above analysis in the algorithmic form 
we have : 
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For J = 1,2,... 

Compute by (4.21) 

Solve for X^ in (a. 2) 

Find the value of Flag, make the correction 
if Flag = 1 until Flag = 0, if X^^^ is close 
enough to X^. 

Compute e(n,s) 

Step 4. Check if u(n,s) satisfied the stopping criteria 
for n = 0,1, ...,N. 

If r* £ r Stop 

If r* > r Update j by j +1 then go to Step 2. 

The application of the above results is seen v;ith the 
illustration in Section V. 
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APPENDIX B 



TABLES 



The contents of this appendix are for the purpose of 
illustrating the convergence of P(n,s), for n = 0,1,..., 6 
as in the example presented in Section V. This is also 
known as a direct application of the results developed in 
Section IV and in Appendix A. 

The values of P(n,s), u(n,s), e(n,s), X^, f(X^) and Flag 
at the j th iteration (j=2,3,***) are tabulated in tables 
numbered from IX to XVIII. 

The following abbreviations are used. 

P(X=n): P(n,s), n = 0,1,.. .,6. 

BOUND: The bound of the error, i.e., u(n,s) 

EST. ERROR: The estimated error, i.e., e(n,s). 
XO: 

FXO: 

IFLAG: Flag. 



The value of X^ . 
f(X^) . 



It is noted that the accuracy desired is achieved at 11 

iterations, with r = .156 10 , the tolerance error in X^ 

-15 -15 

and in f(X^) less than 10 (4)=6 = 10 ), and the maximum 

number of Newton iterations allowed (L*) at 15. 

.The generated tables shown in this appendix and all compu- 
ter work presented in this thesis were done utilizing the 



6l 



grff'-p fe., ^ ' ■ ® ■ wp!i:; : : 

"■^ k-- A. -. *li&*'''' = 

, ^ ■'%*% 










IBM 360/67 computer at the Computer Facility of The Naval 
Postgraduate School, and programmed in Fortran-Formac and 
in Fortran G using double precision. 
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TABLE 



IX 



ITERATION H 2 



PR (X=0) ,eCUND 


0 .375552340 


00 


0.230137360 00 




EST. ERROR 






-C.51C956330-01 




XC,fXC,IFLAG 


0 .112560320 


Ci 


0.138777880-16 


0 


PR (X=n , ECUNC 


0 .378992850 


CO 


0.245171900 OC 




EST. ERROR 






-0.121100290 00 




XC ,FXC ,IFLAG 


0.147353130 


01 


0.0 


0 


PR(X=2) ,eCUNC 


0.251303150 


CC 


C. 176471100 OC 




EST. ERROR 






-C. 110763490 00 




XC, FXC , IFLAG 


0.177509920 


01 


0. 13E777880-16 


0 


PR(X=3), BOUND 


0.134236540 


CC 


0.103553270 OC 




EST. ERROR 






-C. 703672410-01 




XC, FXC, IFLAG 


0.203899420 


01 


0.136777880-16 


0 


PR(X=A) , BOUND 


0 .619228280- 


-Cl 


0.527815110-01 




EST. ERROR 






-0.36 1290140-0 1 




XO,FXO, IFLAG 


0 .227186050 


01 


0.138777880-16 


■ 0 


PR(X=5) , BOUND 


0.255089300- 


-01 


0 .240979980-01 




EST. ERROR 






-C.1598C4C90-01 




XC , FXC , IFLAG 


0 .247886690 


01 


0.190819580-16 


0 


PR(X=6), ECUNC 


C.95677C15D- 


-02 


C. 100345120-01 




EST. ERROR 






-0.629667400-02 




XC, FXC , IFLAG 


0 .266409430 


01 


C.2949C2990-16 


0 


R=» I S 


0.2451719C30 


OC 





GREATER THAN R = 0 . 1 5600 COOOD-05 

ACCLRACY CESIREC HAS NOT BEEN ACHIEVED, 
ITERATION CONTINUES.. 
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TABLE 



X 



ITERATION ft 3 



PR(X=0) ,BCUNC 


0 .416723150 


OC 


0.9903C114D-01 




EST. ERROR 






-0.713398240-01 




XC ,FXG ,IFLAG 


0.10631004D 


Cl 


0.40766CC20-16 


0 


PR(X=l),eCUNC 


0.36744446D 


CO 


0.105499610 OC 




EST. ERRCB 






-0.882C2120D-01 




XC.FXCiIFLAG 


0 .139399430 


Cl 


0.529C9C66D-16 


0 


FR(X=2),ECLNC 


0 .212399260 


00 


0.75937C54D-C1 




EST. ERROR 






-0.5805C477D-01 




XC,FXC tIFLAG 


0.168137400 


01 


C.503C6961D-16 


0 


PR(X=3) ,eCUNC 


0.984821400- 


•01 


0.445598770-01 




EST. ERROR 






-0.27475991D-G1 




XC.FXO , IFLAC- 


0.193306640 


Cl 


-C. 138777880-16 


c 


PB(X=4) ,eCLNC 


0 .391695980- 


01 


0.227123450-01 




EST. ERRCB 






-C. 10 18C979D-01 




XC,FXC,IFLAG 


0.215509190 


01 


0.138777880-16 


0 


FR(X=5),ECUNC 


0. 13 7758050- 


•Cl 


0.103695790-01 




EST. ERROR 






-C.296CC250D-02 




XC,FXC ,IFLAG 


0.235215020 


01 


0.277555760-16 


0 


PB(X=6),ECUNC 


0 .434904960- 


■C2 


0.431793790-02 




EST. ERROR 






-C. 600975350-03 




XC,FXC,IFLAG 


0.252795640 


01 


o 

• 

o 


c 


B=» I S 


0.1054996090 


00 





GREATER THAN B = 0 . 1560CC0CCD-05 

BEEN ACHIEVED, 
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ACCURACY CESIBEC HAS NOT 
ITERATION CONTINUES.. 



TABLE 



XI 



ITERATION # 4 



PR(X=CJ,ECUNC 


0.39622318D 


CC 


0.337685920-01 




EST. ERROR 






-0.296 179920-01 




XCiFXO , IFLAG 


0.10C34698D 


Cl 


C. 303576610- 16 


C 


FR (X=l) , eCUNC 


0.324897360 


00 


0.359746450-01 




EST. ERROR 






-0.233148870-01 




XCtFXO, IFLAG 


0.13170308D 


Cl 


C. 867361740-17 


0 


PR(X=2 ) ,BCUNC 


0 . I7694467D 


OC 


0.25894C16D-01 




EST. ERROR 






-C.77 IC5C91D-02 




XCtFXO, IFLAG 


0.15893C560 


Cl 


C. 277555760-16 


0 


PR(X=3), ECUNC 


0.78528452D- 


Cl 


0.151946130-01 




EST. ERROR 






0.3997C397D-03 




XC,FXC ,If LAG 


0. 18272949D 


Cl 


C. 971445150- 16 


0 


PR (X = 4) , ECUNC 


0.30485588D- 


01 


0.774475420-02 




EST. ERROR 






0.215730340-02 




XC ,FXC , I FLAG 


0 .2036368ID 


01 


0.286229370-16 


C 


PR(X=5), ECUNC 


0.107257C9D- 


•Cl 


0.35359554D-C2 




EST. ERROR 






0.160976930-02 




XCfFXC, IFLAG 


0 .222067690 


01 


0.398986400-16 


0 


PR(X=6), ECUNC 


0.349526330- 


■C2 


C. 147238730-02 




EST. ERROR 






0.835562780-03 




XC.FXC , IFLAG 


0.238344020 


Cl 


0.199493200-16 


0 



R* IS 

GREATER THAN R 



0.359746452D-01 

0.156CCC0CCD-05 



ACCURACY CESIREC HAS NOT BEEN ACHIEVED, 
ITERATION CONTINUES.. 
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8 



TABLE 



XII 



ITERATION U 5 



PR (X=C ) , ECLNC 


0.3B5396460 


00 


0.948962770-02 




EST. EFRCR 






-0.57946C930-02 




XC,FXC ,IFlAG 


0.968632340 


CO 


-C. 520417040-17 


0 


PR (X=l ) , ECLNC 


0 .3.13647170 


00 


0.101095710-01 




EST. ERROR 






-0. 159669190-02 




XC,FXC tlFLAG 


0.127831480 


01 


0.199493200-16 


0 


PR (X=2 ) , ECUNC 


0. 172343960 


CO 


0.727671960-02 




EST. EPRCR 






0.331611240-02 




XC,FXO,IFLAG 


0.155203620 


Cl 


0.43368C870-16 


0 


PR(X=3) , ECUNC 


0.783557160- 


Cl 


0.426998C30-02 




EST. ERROR 






0.301330580-02 




XC ,FXC ,IFLAG 


0.179683850 


01 


0.312250230-16 


c 


PR{X=4), ECUNC 


0. 315700990- 


Cl 


0.217642580-02 




EST. EPRCR 






0.166560150-02 




XC ,FXC , IFLAG 


0 .201837060 


0 1 


C.549 148400- 16 


0 


PR(X=5), ECUNC 


0.116401170- 


Cl 


0.993671890-03 




EST. EFRCR 






0.688320030-03 




XC ,FXC , IFLAG 


0.222124790 


Cl 


0.286229370-16 


0 


PR (X=6) , ECLNC 


0 .399688130- 


■C2 


0.413769320-03 




EST. EFRCR 






0.223634820-03 




XC.FXO , IFLAG 


0.240924100 


01 


0.294902990-16 


c 


R=» IS 


0. 101C957 140- 


-0 1 





GREATER THAN R = 0. 156CCC0000-05 

ACCURACY CESIREC HAS NOT BEEN ACHIEVEDi 
ITERATION CONTINUES.. 
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I 



ITERATION # 



6 



PR(X = C) ,eCLNC 
EST. ERROR 
>C,FXC ,IFLAG 
PR (X=l ) , ECUNC 
EST. ERROR 
XCfFXG .IFLAC- 
PR (X = 2 ) ,BCUNC 
EST. ERROR 
XO.FXO ,IFLAG 
PR (X=3 ) , EOLKC 
EST. ERROR 
XC ,FXO , IFLAG 
PR (X = A ) , BOLNC 
EST. ERROR 
XO.FXC , IFLAG 
PR(X = 5 ) .EGLNC 
EST. ERROR 
XO,FXO , IFLAG 
PR(X=6), BOLNC 
EST. ERROR 
XO,FXC , IFLAG 

R=fr IS 

GREATER THAN R 



0.383426950 00 

C.95860549C CC 
0.312240380 OC 

0.127401480 01 
0.173358040 CO 

0.155709650 01 
0.794206330-01 

0.181331160 01 
0.321835230-01 

0.204653790 01 
0.118926590-01 

0.225948410 01 
0.407474920-02 

0.245407260 01 



0. 226213660-02 
0.251575800-06 
0.199493200-16 
C.24C9S2C80-02 
0.159803520-02 
-0.303034510-16 
0. 172462530-02 
0.141573780-02 
0.525838050-16 
0.101787840-02 
0.618948720-03 
0.468275340-16 
0.518816630-03 
0.129078450-03 
C. 867361740-18 
0.236871620-03 
-0.267 17 1500-04 
-0.329597460-16 
0.986342780-04 
-0.392665630-04 
C. 149186220-15 



0. 2409920810-02 
0.1560000000-05 

BEEN ACHIEVEC, 



ACCLRAGY CESIREC HAS NOT 



ITERATION CONTINUES.. 



TABLE 



XIV 



ITERATION # 7 



PR (X = C) ,ECLN'C 


0 .383289720 


CO 


0.46752E92D-03 




EST. EPRCR 






0.28859S40D-03 




XC.FXCtIFLAG 


0.958211450 


00 


-0.384891770-17 


0 


PR (X=l ) , ECUNC 


0.3^3491610 


CC 


0.498C7191D-03 




EST. ERROR 






0.402426830-03 




XC,FXC tIFLAG 


0.127701890 


Cl 


0.356328820-16 


0 


PR(X=2) , ECUNC 


0. 173648390 


CO 


0.3585C477D-03 




EST. ERROR 






0.124923690-03 




XC ,FXO tIFLAG 


0.156283720 


Cl 


C. 254245410-16 


0 


PR (X = 3) , ECUNC 


0.795506890- 


Cl 


0.21C37C66D-03 




EST. ERROR 






-C. 457962140-04 




XC.FXC tlFLAG 


0.182001440 


01 


C. 108420220-1 7 


0 


FR(X=4 ), ECUNC 


0.322014900- 


Cl 


0. 107226760-03 




EST. ERROR 






-0.644269500-04 




XC ,FXC ,i FLAG 


0.205198530 


Cl 


0.241234980-16 


0 


FR(X = 5 ), ECUNC 


0.118764480- 


01 


C. 489555920-04 




EST. ERROR 






-C. 360735770-04 




XC ,FXC .IFLAG 


0.226164110 


01 


C.99628C15D-16 


0 


PR(X = 6) , ECUNC 


0 .^05946180-02 


C. 203653230-04 




EST. ERROR 






-0.137C9141D-04 




XC,FXC , IFLAG 


0.245153090 


01 


0.33556C57D-15 


C 


IS 


0.498C71909D 


-03 





GREATER THAN R = 0. 156C0COCCD-O5 

ACCURACY CESIREC HAS NOT BEEN ACHIEVED, 
ITERATION CONTINUES. . 
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TABLE 



XV 





ITERATION fli 6 








PR(X=OJ,eCUNC 


0 .383296490 


CC 


0.852C87360-04 




EST. ERROR 






0.724810970-04 




XC ,FXC , IFLAG 


0.958757370 


00 


C. 454348470- 16 


C 


PR(X=1 ) ,ECUNC 


0.3.13535800 


CO 


C. 907753000-04 




EST. ERROR 






0.259713920-04 




XC,FXC, IFLAG 


0.127807790 


Cl 


C.954C97910-16 


0 


PR(X=2) , BOUND 


0 .173664180 


CC 


0.653387150-04 




EST. ERROR 






-0.304558580-04 




XC.FXC, IFLAG 


0.156364080 


01 


-C. 461124740-16 


0 


PR (X=3 ) ,ECLNC 


0.795414800- 


Cl 


0.3834C7690-04 




EST. ERROR 






-0.2983C7340-C4 




XC,FXC, IFLAG 


0.181992940 


Cl 


C. 310152970-15 


0 


PR(X=4i .BOUND 


0.32 1895020- 


■C L 


0.195424410-04 




EST. ERROR 






-0. 1263362C0-04 




XC.FXC tlFLAG 


0 .205074740 


01 


0.381256310-15 


C 


PR(X=5), BOUND 


0 . 1.18697910- 


Cl 


0.892232310-05 




EST. ERROR 






-0.252257130-05 




XC ,FXC , IFLAG 


0.225937230 


Cl 


0.147939390-15 


0 


PR (X=6i , BOUND 


0 .405705590- 


C2 


0.37 1529440-05 




EST. ERROR 






0.430411580-06 




XC.FXC .IFLAG 


0 .244863480 


01 


C. 639679280- 17 


C 



R* IS 0. 9077529970-04 

GREATER THAN R = 0 . 1560CC0C0D-05 

/CCURACY CESiREC HAS NOT BEEN ACHIEVED, 
ITERATION CONTINUES. . 



69 



TABLE 



XVI 



ITERATIOK # 9 



FR (X = C ) , ECUNC 


0.38329881D 


CO 


0.138EC549D-04 




EST. EFRCR 






C.772C1441D-05 




XC,FXC ,IFLAG 


C.95892C26D 


CC 


0.85665E530- 17 


C 


FR (X=l ) ,ECLNC 


0.31353855D 


cc 


0.147E7346D-04 




EST. ERRCR 






-C.68844768D-05 




XC,FXC ,IFLAG 


0 .12781740D 


Cl 


0.728E8550D-16 


0 


FR(X=2),ECLNC 


0 .17366227D 


CO 


0.10643712D-04 




EST. ERRCR 






-0.85494961D-05 




XC,FXC tIFLAG 


0.15635257D 


Cl 


C.11780534D-16 


0 


FR(X = 2),eCL'NC 


0.795390UD- 


Cl 


0.62457321D-05 




EST. ERRCR 






-0.26 186557D-05 




XC,FXC ,IFLAG 


0.18196I82D 


Cl 


C.27816562D-16 


0 


PR (X=4 ) , ECLNC 


0.32168494D- 


Cl 


0.31834742D-05 




EST. EFRCR 






C.53571488D-06 




XC.FXC ,IFLAG 


0.20503695D 


01 


-0.81416807D-17 


0 


PR (X = 5 ) , ECUNC 


0 . 11869734D- 


Cl 


0. 14534513D-05 




EST. EFRCR 






0.85567C10D-06 




XC,FXC flFLAG 


0 .22590925D 


Cl 


C.762l^280D-l6 


c 


PRIX=6), ECLNC 


0.40572299D- 


02 


0.60522347D-06 




EST. EFRCR 






0.43849225D-C6 




XC,FXC,IFLAG 


0 .24485899D 


Cl 


-0. 162996C6D-15 


c 


R* I S 


0.i47873455D- 


-04 





GREATER THAN R = 0 . 1 56 00C0 CCD-0 5 

/CCLRACY CESIREC HAS NOT BEEN ACHIEVED, 
IlERATICN CCNTINLES. . 
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TABLE 



XVII 



ITERATION U 10 



FR(X=0) ,eCLN'C 


0.383299C60 


CO 


0.204332300-05 




EST. EFRCR 






-C.1591C869D-06 




XC,FXG,IFLAG 


0 .95894C18D 


CC 


-0.243675610-16 


0 


FR(X = n .ECUNC 


0 .313538530 


CO 


0. 217681040-05 
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